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We compare three techniques to extract partial photoelectron spectra from the wave packet re- 
sulting from the integration on a finite-element discrete variable representation basis of the time- 
dependent Schrodinger equation for an helium atom subject to an ultra-short XUV pulse. These 
techniques are: projection on products of hydrogenic bound and continuum states, projection onto 
multi-channel scattering states computed in a B-spline close-coupling basis, and a technique based 
on exterior complex scaling (ECS) [Palacios et al, Phys. Rev. A 76, 043420 (2007)] implemented 
in the same basis used for the time propagation. These methods allow to monitor the population 
of continuum states in wave packets created with ultrashort pulses in different regimes. The first 
method works well at the energies where the ionization continuum is unstructured while it becomes 
inefficient close to threshold openings due to the presence of long-hved metastable states or of van- 
ishingly slow photofragments. The agreement between the last two methods is excellent below the 
double ionization threshold. This result shows that in regions where the projection on Coulomb 
functions is inapplicable, scattering states computed in an optimized basis can be converted to a 
representation well suited to time-dependent simulations. Above the double ionization threshold, 
the projection onto Coulomb functions and the ECS method are found in excellent agreement with 
each other. 

PACS numbers: 31.15.ac, 32.80.Fb, 32.80.Rm, 32.80.Zb 



I. INTRODUCTION 

During the last decade, two transformational experi- 
mental techniques, high harmonic generation [1] and x- 
ray free electron lasers [2] , have given access to femtosec- 
ond and sub-femtosecond intense light pulses in the XUV 
and soft x-ray energy range, thus opening the way to time 
resolved studies of the correlated motion of electrons in 
atoms and molecules on their characteristic time-scale [3- 
10]. These new techniques can be used not only to mon- 
itor the electronic motion but also to steer it [11, 12]. 
This latter capability offers the perspective of control- 
ling dynamics at the femtosecond and sub-femtosecond 
timescale, such as electronic dynamics in atoms [13, 14], 
molecules [15-17], and solids [18], and eventually also nu- 
clear dynamics such as fast proton migration [19]. 

The interpretation of experiments on attosecond dy- 
namics, however, face a number of difficulties and re- 
quire guidance by theory. First, sub- femtosecond pulses 
typically excite the target to a coherent superposition of 
states above the ionization threshold and across a wide 
range of energies. As a consequence, several different 
ionization regimes such as multiply excited autoionizing 
states, multichannel single ionization states and, possi- 
bly, multiple ionization states, are accessed at the same 
time. Second, in common pump-probe schemes [15], the 
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strong few-cycle IR-probe pulse that follows an attosec- 
ond weak XUV-pump pulse gives rise to electronic dy- 
namics that unfolds on a very short time-scale through 
non-perturbative stages, e.g. tunneling, over-the-barrier 
ionization, and Rabi oscillations. 

Traditional perturbative approaches [20] are clearly 
not well suited to describe such dynamical regimes. 
Moreover, since the duration of compressed IR pulses 
easily spans just a few [21] or even a single [22] carrier 
cycle, stationary non-perturbative techniques like those 
based on the Floquet method [23] cannot be used either. 
Reliable theoretical predictions for ultrashort processes, 
therefore, generally require direct integration of the time- 
dependent Schrodinger equation (TDSE) [14, 24-34]. 
Such an approach permits to reproduce faithfully the 
physical process under study. However, it gives rise 
to a multitude of problems as well. Relevant parts, if 
not most, of the electron dynamics triggered by sub- 
femtosecond pulses take place in the ionization contin- 
uum. Indeed, typical experiments are designed to mon- 
itor the energy and angular distribution of the photo- 
electrons emerging from the reaction center, e.g., with 
a velocity map imaging spectrometer [35] or even sev- 
eral photo- fragments in coincidence, e.g., with a reaction 
microscope [3G]. Transient absorption spectroscopy [37], 
which monitors quasi-bound electronic dynamics, consti- 
tutes a notable exception. One of the most prominent 
problems theory has to face is then how to extract from 
a numerical simulation, intrinsically limited in both time 
and space, the relevant asymptotic scattering informa- 
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FIG. 1. Approximate regions of applicability of different 
extraction methods for the single ionization continuum of he- 
lium as a function of photon energy and pulse duration for 
a single-photon transition. SI: single ionization, DE: double 
excitation, DI: double ionization. Green shaded: projection 
on scattering states, blue solid: projection on products of 
Coulomb functions, red solid: Berkeley ECS method. See 
text for details. The upper part shows the (logarithmic) ion- 
ization probability for single ionization and positions of ex- 
cited bound states. 



tion. 

A number of alternative techniques have been used 
in the past to extract differential distributions of the 
products of a light-induced reaction from an entan- 
gled wave function: spectral analysis of autocorrelation 
functions [38], analysis of the radial flux at large dis- 
tances [39], projection of the wave packet onto products 
of one-particle functions, e.g., Coulomb functions ['2G, 
32, 40] or Volkov states [28], surface integration of 
the monochromatic component of the wave function ex- 
tracted with a technique based on exterior complex scal- 
ing (Berkeley-ECS) [41, 42], or projection on scattering 
states computed on the same basis used to carry out the 
time propagation [14, 29-31, 40, 43, 44]. 

In the following, we compare three of these techniques 
in more detail for the prototypical three-body system, the 
helium atom: projection on products of one-particle func- 
tions, projection on scattering states, and the Berkeley- 
ECS method. For all these techniques, the extraction 
of asymptotic scattering information takes place in the 
field-free region after the laser pulse is over. As will be 
discussed in more detail later, each of these approaches 
has its strengths and drawbacks making it applicable for 
different photon energies and pulse durations. These re- 
gions are schematically illustrated in Fig. 1. 

We here give a brief preview. The projection on 
products of Coulomb functions works well when the 
wave packet only contains electrons that are already 



well-separated in radial space. While this is relatively 
straightforward to achieve at high energies (such as above 
the double ionization continuum) , it becomes prohibitive 
close to the doubly excited resonances that populate the 
energy region from around —0.7 a.u. to Oa.u., and which 
only decay after many femtoseconds. The main advan- 
tages of this approach are its analytical and numerical 
simplicity and the fact that it works above the double 
ionization threshold. This allows it to be used even when 
very large computational boxes are employed, and con- 
sequently with very long driving pulses. In the case of 
single ionization, a closely related approach with almost 
identical properties is the projection on channel func- 
tions, i.e., eigenstates obtained by freezing the inner elec- 
tron in an ionic state. 

Projection on scattering states (PSS), on the other 
hand, works for any separation of the final particles, and 
thus can be used also when doubly excited resonances are 
excited and before they have decayed. The algorithms 
to calculate atomic and molecular scattering states have 
now achieved a considerable level of sophistication and 
can be used with large radial boxes, i.e., for large sep- 
aration of the electrons. Moreover, the calculation of 
scattering states is independent of the actual simulation 
and can be optimized separately. Yet, scattering states 
have so far been employed only in those cases where they 
are built in the same basis which is used to carry out the 
simulations. We demonstrate below that the scattering 
states computed in an optimized B-spline close-coupling 
basis with the K-matrix method [45] can be accurately 
converted to a finite-element discrete-variable (FEDVR) 
basis [32] optimized for time propagation. Fully differ- 
ential photoelectron spectra can thus be extracted from 
the wave packet arising in photoionization by means of 
a simple projection. However, scattering states are not 
available above the double ionization threshold as the 
boundary conditions for double ionization [46] cannot be 
easily enforced due to the infinite set of constraints they 
entail (as opposed to the finite number of constraints of 
single ionization problems). In addition, the scattering 
states become exceedingly expensive to calculate as the 
double ionization threshold is approached from below due 
to the presence of many double Rydberg series. 

Finally, the Berkeley ECS method is applicable for 
both single and double ionization, and also works when 
resonances and other long-lived states are involved. Its 
main drawback is the computational complexity: for each 
desired final energy, a large linear system representing the 
exterior-complex-scaled Hamiltonian acting on the final 
wave packet has to be iteratively solved. This approach 
becomes computationally expensive when the system has 
large spatial extent and when the wave packet to be an- 
alyzed covers a wide range of energies. With currently 
available supercomputers, this becomes impractical for 
linear systems with a dimension of more than a few mil- 
lion. For typical simulations in helium this limits the 
box sizes to the range of ~ r^ax = 250 a.u. and thus to 
simulations where only relatively short pulses are used. 
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For the comparison of the methods, we choose the 
pulse parameters to highlight the limitations of the 
"most straightforward" projection methods, the projec- 
tion on the product of Coulomb functions or similar 
single- channel functions. These limitations emerge most 
markedly when autoionizing states are excited, in which 
case convergence cannot be achieved. Such behavior is 
to be contrasted with the Berkeley-ECS method and the 
projection on scattering states which yield convergent re- 
sults in excellent agreement with each other. We demon- 
strate that scattering states can be fruitfully ported 
across different bases, thus giving access to a numeri- 
cally convenient and accurate method of monitoring the 
single ionization component of a wave packet below the 
double ionization threshold. 

The article is organized as follows. In section II we 
introduce the different approaches that are used in this 
work to extract physical information from a wave packet 
and put them in comparison from a theoretical perspec- 
tive. In section III we briefly summarize the propaga- 
tion method for the time dependent wave function in the 
finite-element discrete- variable representation. In section 
IV we illustrate the details of the methods for extracting 
information from the time-dependent wave packet. In 
section V we present simulation results highlighting the 
similarities and differences between different extraction 
methods. Finally, in section VI we draw our conclusions. 



II. THEORETICAL PRELIMINARIES 

The electronic dynamics induced in an atom or 
molecule by attosecond pulses is encoded in the energy 
and angular distribution of the emergent photofragments. 
To unravel the mechanism underlying the break-up pro- 
cess it is thus necessary to establish a correspondence be- 
tween these experimental observables and the simulated 
wave packet dynamics. In principle, this is straightfor- 
ward: since the photofragments are collected at macro- 
scopic distances from the reaction center, the final kine- 
matic state of the system is fully identified by a com- 
plete set of quantum numbers for the separated frag- 
ments. In practice, however, several complications arise 
which we briefly review in the following. Upon conclu- 
sion of the electromagnetic pulse, the numerical solution 
of the TDSE yields the wave function '^/(t) of the wave 
packet. The dynamics of \E'(i) is then governed by the 
total field- free atomic Hamiltonian, Ha- The state ^'(t) 
is composed of a bound and an unbound part. Let us 
indicate by A the projector on the bound states of Ha, 
and with ^''(t) = (1 - A)4'(t) the unbound part of ^'(i). 
The fragments detected are associated with the long-time 
limit of the unbound component ^'{t). It consists of 
a superposition of unbound eigenstates ipae of the sum 
of the Hamiltonians of the separated fragments, Hq. In 
mathematical terms the latter is frequently referred to as 



asymptotic completeness [47, 48] 
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In Eqs. (1,2), a designates the collective set of quantum 
numbers beyond the total energy identifying uniquely the 
final state of the system. These include the fragmentation 
channel of the target, the asymptotic angular distribu- 
tion of the photofragments, and their internal quantum 
numbers (e.g. spin and angular momentum). 

The expansion coefficients Caeit) in (1) may not nec- 
essarily converge in the infinite-time limit. In particular, 
when two or more fragments are charged, the phase of 
Cae{t) diverges due to the long-range character of the 
Coulomb field. Notwithstanding this difficulty the prob- 
ability density |cQe(i)p will still converge in the sense 
of distributions (i.e., when convoluted with a test func- 
tion) yielding a well defined asymptotic distribution of 
the fragments. 

In the case of two charged fragments, a stronger result 
is known to hold: the absolute value |cQ,e(t)| itself con- 
verges [48] 'AS t oo. Hence, the probability distribution 
Pa{E) for the detection of fragments in channel a and 
with energy E at the end of the propagation is given by 



P„(S) = lim \c^E[tt 



(3) 



More details on this result can be found in Sec. XI.9 
of [48]. Moreover, a substantial body of numerical evi- 
dence [30, 34, 40, 49, 50] indicates that the more lenient 
result (3) holds in the case of atomic double photoioniza- 
tion as well. 

For the time being we will thus assume that, with suit- 
able precautions, Eq. 3 is applicable. The probability 
distribution Pa (E) can then be computed as 



Pa{E) = lim 

t^OO 



(4) 



At first sight, the prescription (4) to extract experimen- 
tal observables from a wave packet has the appeal of sim- 
plicity since the uncoupled states (paE are usually more 
easily obtainable than the continuum eigenstates of the 
full Hamiltonian. This simplicity, though, is misleading 
since the projector A of the total Hamiltonian requires at 
least a certain number of bound states of the fully inter- 
acting system to be known. Since the bound states of Ha 
are not orthogonal to (paE, their elimination is essential. 
Otherwise they would give rise to spurious contributions 
to the ionization channels which do not vanish for large 
times. 

The greatest drawback of Eq. 4, however, stems from 
the asymptotic time limit, which may be out of reach for 
the propagation algorithm. Even by the inclusion of part 
of the long-range interactions between photofragments 
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into the reference Hamiltonian Hq this problem can only 
be marginally alleviated rather than solved. More gener- 
ally, all the methods that require the wave packet to reach 
the asymptotic region, i.e. the region where the dynam- 
ics governed by Hq and Ha become equivalent, face the 
same problem, namely the propagation of the fully corre- 
lated wave function for long times and at large distances 
which may be computationally prohibitively expensive. 
This limitation becomes particularly severe in a num- 
ber of circumstances frequently encountered in atomic 
photoionization. For example, when the wave packet 
spectrum is distributed across an ionization threshold, 
the wave packet comprises components with vanishingly 
small kinetic energy that take exceedingly long times to 
reach the asymptotic region. Even more severe, when 
several channels with different thresholds are simultane- 
ously open, both slow and fast photoelectrons are present 
at the same time. Hence, in order for the slowest part 
of the wave function to reach the asymptotic region, the 
propagation box must be large enough to accommodate 
for the fastest components as well. A further difficulty 
with Eq. 4, perhaps the most relevant in the present con- 
text, is provided by resonances in general and by Rydberg 
series of doubly excited states in particular. Such doubly 
excited states are a general feature in atomic photoioniza- 
tion spectra. First, the convergence of the resonant pro- 
files in those channels where the excited resonances decay 
requires a propagation time proportional to the lifetime 
of the longest lived resonance which is excited in the sim- 
ulation. Second, doubly excited states have, in general, 
non- vanishing scalar products with all the eigenfunctions 
of the unperturbed Hamiltonian, including those belong- 
ing to closed ionization channels. 

One avenue to circumvent some of the limitations as- 
sociated with Eq. 4 is to use, instead of the asymptotic 
states (faE, the scattering states of the total atomic 
Hamiltonian Ha [47]. The (j>~^ states fulfill so-called in- 
coming boundary conditions [51, 52], which are appro- 
priate to the context of photoionization, since photoelec- 
trons are observed in the positive time limit. 

Such scattering states satisfy the Lippmann-Schwinger 
equation with advanced Green's functions [47] 

^-E^Vc.E + G^iE)H'^p-^ (5) 
= ^a.E + G-{E)H'^^E (6) 

where the operators Gq{E) — {E — Hq — iO+)^^ and 
G~{E) = {E—Ha — iO'^)~^ are the resolvents of the refer- 
ence {Hq) and full Hamiltonian {Ha), while H' = Ha—HQ 
is the corresponding perturbation (the interactions not 
accounted for by the channel Hamiltonian Hq). 

The ip~E states form a complete orthonormal basis for 
the unbound states of Ha, 

Ha^p-E^E^~E: i^p-^l^p-^,) ^ S^pS{E - E') (7) 
l-A = J2 [ de\^-,){^-,\. (8) 



We can thus write the scattering component '^'{t) of the 
wave packet at any time t after the external field is over 

as 

\^i>'[t)) = e-'^°(*-*«)(l - A)|*(to)> = 

= E / rfe|CJe-"(*-*°)c^,(io) (9) 

where 

Cai;(^0) = {^-E\^'{tQ)). (10) 

A crucial aspect of equation (9) is that, in the large time 
limit, the states can be replaced by their asymptotes 
ipae [51, 52]. This is one defining feature of the -0" states, 
sometimes referred to as control of 0~ by the future. In 
particular, 

lim |c„,(t)| = |c- (to)| (11) 

r— >C30 

To compute the distribution Pa{E) we now com- 
bine Eq. 3, Eq. 11, and Eq. 10 and obtain the exact result 

Pc.{E) = \{i^-Emto))\'- (12) 

In contrast to Eq. 4, Eq. 12 requires neither a projection 
on the bound states of the system, to which the scatter- 
ing states are orthogonal, nor the evaluation of a long- 
time limit. The convenience of using scattering states 
ipaE^ instead of the asymptotic limits (paE, is thus ap- 
parent: Eq. 12 can be evaluated as soon as the external 
time-dependent field is over, without having to wait until 
the ionizing wave packet reaches the asymptotic region 
(Eq. 4). 

In [40] , the use of scattering states has been deemed im- 
practical but in the simplest cases, due to the additional 
workload required to solve Eq. 5 or Eq. 6. Indeed, no sys- 
tematic procedure to generate accurate scattering states 
for the double ionization of atoms has been reported to 
date. As far as single-ionization processes are concerned, 
however, the calculation of scattering states is straight- 
forward. This task has been tackled successfully in the 
course of the last four decades [53-57]. Today several effi- 
cient methods capable of computing multi-electron single 
ionization scattering states are available. They include 
the R-matrix [58], J-matrix [59], K-matrix [45, 60], Fesh- 
bach projection [01], and inverse iteration [62]. Further- 
more, the computational overhead for generating scatter- 
ing states is easily compensated whenever a large number 
of different simulations must be carried out, as is the case, 
e.g., of time-delay scans in pump- probe schemes. The 
biggest drawback of these methods is that they only work 
for total energies below the double ionization threshold. 

A third elegant strategy for computing Pa{E), based 
on exterior complex scaling, was put forward by Pala- 
cios, McCurdy and Rescigno [33, 41, 63]. In this ap- 
proach, referred to in this work as the Berkeley-ECS 
method, the monochromatic component at energy E 
of the wave packet 4'(io) is extracted by applying to ^'(to) 
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the retarded resolvent G~^{E) of the total Hamiltonian, 
'i/E = {E)'i> (to) , the realization of which in a finite 
radial domain is provided by the resolvent of an exterior- 
complex- rotated Hamiltonian Hg . From the function ^ e 
the distribution of the photofragments is extracted in the 
asymptotic region. This last step is accomplished [33, 63] 
by carrying out a surface integral which involves 5* e and 
test functions analogous to the (paE states introduced at 
the beginning of this section (see Sec. IV). 

The Berkeley-ECS method is accurate [33, G4] and is 
not restricted to single ionization. It can be implemented 
directly in the same basis that is used to carry out the 
integration of the TDSE, irrespective of the number of 
electrons that are eventually produced. With the super- 
computers available today, such an approach is feasible 
for moderate box sizes but prohibitively expensive for 
larger problems. Both the projection on scattering states 
(PSS) and the Berkeley-ECS method overcome some of 
the limitations of methods based on the use of approx- 
imated continuum functions while suffering from some 
drawbacks of their own. 

On the one hand, the Berkeley-ECS method is appli- 
cable at energies above the double ionization threshold. 
On the other hand, the PSS method is more convenient 
below the double ionization threshold. This is because 
the single ionization space is only a tiny fraction of the 
discretized continuum which is apt to represent multiple 
ionization processes. The calculation of scattering states 
takes advantage of this reduced state space by employing 
compact close-coupling representations. 

Moreover, the two methods have different scaling prop- 
erties. While scattering states are computed once and for 
all, and can thus be applied to the analysis of an arbi- 
trary number of wave packets through the evaluation of 
a simple scalar product, the time-consuming steps of the 
Berkeley-ECS method must be repeated each time a wave 
packet is to be analyzed. 

The use of scattering states furthermore permits to 
monitor and isolate the different channel components of 
the wave functions during the time evolution, in par- 
ticular to disentangle individual ionization pathways in 
pump-probe experiments. 

The latter holds, of course, only if scattering states 
can be converted from a basis best suited to their effi- 
cient calculation to another one which is best suited to 
integrate the TDSE. We will show in the following that 
this is indeed possible. 



III. PROPAGATION METHOD 

Our computational approach (see [32, 65, 66] for a 
more detailed description) for solving the time-dependent 
Schrodinger equation for two-electron systems, 

i|-*(ri,r2,t) = i/*(ri,r2,t), (13) 



is based on a time-dependent close-coupling scheme 
[67-70] where we expand the angular part of the six- 
dimensional wave function *(ri,r2) in coupled spherical 
harmonics J;^^;^ (f2i, ^la)- 

The interaction of a helium atom with linearly polar- 
ized light is described by the Hamiltonian 

Tir _ rr I ttI.v _ Pi P2 2 2 1 , 

if - + --+----- -t- ^^—^ + , 

(14) 

where the interaction with the electromagnetic field in 
the dipole approximation, H^J^, is either given in length 
or velocity gauge. The gauge independence of the phys- 
ical observables is a necessary condition for the conver- 
gence of the numerical solution. 

For the discretization of the radial functions 
Rf" I (ri ,r2,t), we employ a finite-element discrete- 
variable representation (FEDVR) [71-73]. We divide the 
radial coordinates into finite elements in each of which 
the functions are represented in a local DVR basis 

with a corresponding Gauss-Lobatto quadrature to en- 
sure the continuity of the wave function at the element 
boundaries. This method leads to sparse matrix repre- 
sentations of the differential operators and to a diagonal 
potential matrix (within quadrature accuracy), enabling 
efficient parallelization. 

For the temporal propagation, we employ the short it- 
erative Lanczos (SIL) method [74-76] with adaptive time- 
step control. The initial He ground state l^S(ls^) is ob- 
tained by relaxing an arbitrary test function in imag- 
inary time. For an initial 2-'^S(ls2s) metastable state 
we directly solve the eigenvalue problem of the field-free 
Hamiltonian (14) in a small box using the SLEPc library 
[77]. 



IV. EXTRACTION OF THE SPECTRUM 

The four methods used here to extract the single ion- 
ization photoclectron spectrum arc: (i) projection of the 
wave packet on the product of a bound hydrogenic wave 
function with Z=2 and a Coulomb continuum function 
with Z=l; (ii) projection on channel functions obtained 
from the diagonalization of the total Hamiltonian in the 
configuration space where one of the two electrons is 
frozen in an hydrogenic state of the parent ion. Channel 
functions differ from the products of Coulomb functions 
because, in the radial region where the bound electron 
density is not zero, the effective potential felt by the free 
electron departs from that of a nuclear charge with Z=l. 
Yet, at larger distances, these channel functions converge 
to phase-shifted Coulomb functions. Therefore, for the 
purpose of evaluating the absolute value of the projec- 
tion of an outgoing wave packet, channel functions and 
pure Coulomb functions are essentially equivalent; (iii) 
projection on scattering states (PSS) obtained with the 
B-spline K-matrix method; and (iv) extraction of the 
partial differential photoclectron distribution with the 
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Berkeley-ECS method. For double ionization we com- 
pare the Berkeley-ECS method with the projection of the 
wave packet on a symmetrized product of two Coulomb 
functions with Tj—1. 

Our implementation of the Berkeley ECS method 
closely follows that of Palacios et al. [33, 63], which builds 
on earlier work reviewed in [78]. We give a brief review 
here and refer the reader to the original papers for de- 
tails. The central idea is to solve the inhomogeneous 
linear system 

{E~Ha)\-^sciE)) = \^{ta)) (15) 

in the basis used for temporal propagation, using the 
PETSc package [79]. The scattered wave '^sc{E) 
is equivalent to the Fourier transform of the time- 
dependent wave packet from t = to to t = oo, and 
corresponds to the application of the retarded Green's 
function of the total Hamiltonian on the wave packet, 

|M'sc(i?)) =G+(i?)|vI'(io)). (16) 

Purely outgoing boundary conditions are enforced by an 
exterior complex scaling (ECS) transformation for each 
of the radial coordinates. As the wave packet at the end 
of the pulse is a square-integrable function, the asymp- 
totic behavior of the scattered wave can be deduced from 
the asymptotic form of the Green's function. For single 
ionization, the amplitude CaE{to) can be expressed as 
[78], 

CaEih) = {Vc.e\E - Ho\^sc{E)}, (17) 

where ipaE is an eigenstate of the reference Hamiltonian 
Hq. The latter should contain the monopolar long-range 
interaction between the fragments to suppress spurious 
contributions. Using Green's theorem allows one to ex- 
press the single ionization amplitude as a surface integral 
in the non-scaled region of space, 

cMto) ^ \ [ ((^aBV*sc(^^) - *Sc(S)V</J„b) • dS, 
•J s 

(18) 

where V = (Vi, V2) is the six-dimensional gradient op- 
erator (in the present case of helium). Since the inte- 
gral (Eq. 18) is evaluated in a radial region far from the 
atom and ^gc(£') is an outgoing wave packet by construc- 
tion, it is sufhcient that (j)aE satisfies the same outgoing 
boundary conditions as the eigenstates of Hq- Thus, in 
single ionization, ipaE can be taken as the symmetrized 
product of an ionic bound state and a Coulomb wave 
with Z = 1. For double ionization, a similar expression 
can be found [G3, 78]. 

The single-ionization scattering states of helium be- 
low the double ionization threshold of the atom are com- 
puted with the B-spline K-matrix method. B-splines are 
a convenient tool to accurately represent the radial com- 
ponent of continuum atomic orbitals on a finite inter- 
val [80, 81], while the K-matrix method [82] is an 
realization of configuration interaction in the continuum 



along the lines of Fano's pioneering paper [83]. The 
K-matrix method has been successfully applied for the 
single-photoionization spectrum of several atomic and 
molecular systems [45, 56, 84-86]. We will provide here 
only a brief description of its implementation for the case 
of helium (details can be found elsewhere [45, 87]). 

A complete set 4'aE stationary eigenfunctions of the 
field-free Hamiltonian Ha at a given energy E in the 
single-ionization continuum are sought in the form of a 
linear combination of partial-wave channel functions c/jaE 
(PWC's), plus an additional component from a localized 
(or pseudostate) channel (LC), 

i^aE = ^l^aE +^^d€(j)^^-I^^K^^^aE , (19) 

7 

where the index a runs over the channels which are open 
at energy E, while the index 7 runs over all open and 
closed channels, including the localized one. 

The PWC a is defined by coupling and antisymmetriz- 
ing a bound state of the He"*" parent ion with quantum 
numbers Na and and energy Ea , to an electron state 
with angular momentum la , the radial degree of freedom 
of which is otherwise unconstrained, to give a state with 
definite values for the total spin S and angular momen- 
tum L. 

<Pc.E = AQs^yltl^^i,^2)RN^LAn)^-^^^, (20) 

where A is the antisymmetrizer, 0ss is a two-electron 
spin function, RN^La is the radial part of the frozen He'^ 
parent ion state, and f^E the continuum radial func- 
tion. Asymptotically the faE are a linear combination 
of the regular and irregular Coulomb functions with an- 
gular momentum ^q, energy E, and a phase shift 5aEi 
determined by the short range behavior of the differen- 
tial equation for faE, which differs from that of the hy- 
drogenic functions. This results from the deviation of 
the frozen-core potential from that of a pure Coulomb 
potential. 

The PWCs in (Ec^. 19) do not exhaust the state space 
associated to single ionization, because the set of bound 
states of the parent ion is not complete, and because the 
close-coupling expansion (Eq. 19) is truncated. Neverthe- 
less, if all single- and double-ionization closed channels 
were to be included in the close-coupling expansion, their 
contribution would decay exponentially at large radii. 
Therefore, instead of using a complete basis, it is suf- 
ficient to include in (Eq. 19) a pseudo-state channel LC 
that comprises a sufficiently large number of normalized 
two-electron functions built from localized orbitals, to 
attain good accuracy. 

Equation (19) may be solved for the unknown coeffi- 
cient matrix K by requiring il)^;^ to be an eigenfunction 
of the complete projected Hamiltonian with eigenvalue 
E, 

{<l>pE-\E~Ha\iZE)^^ ^fi,E'. (21) 



This condition leads to a system of integral equations for 
K which can be discretized and solved with standard lin- 
ear algebra routines. The scattering states with definite 
spherical symmetry are then computed as 

(22) 

where 'K.afiiE) = J^aE./SE is the on-shell reactance ma- 
trix (§7.2.3 in [47]) while ae^ and da are the Coulomb 
and channel phase shifts, respectively. Finally, the scat- 
tering states which correspond to Coulomb plane waves 
associated to a parent ion in a given state A, are given 
by 

JVo,=N_4 

a 

where L^, Mj^, and indicate the parent-ion's angu- 
lar momentum and spin, f2 and a indicate the asymp- 
totic photoelectron's direction and spin, and C^^ are 
Clebsch-Gordan coefficients. The states ^/i^ are nor- 
malized according to 

ii^lEnJ^B.E'n'.') ^ SAB5aa'S{E- E')S{n-n'). (24) 



V. NUMERICAL EXAMPLES 

In this section we present results of alternative extrac- 
tion methods applied to helium under the action of an 
ultrashort XUV laser pulse, 

Re + ^ lle+* + e- . (25) 

Helium is the simplest system that features autoionizing 
states and excited-threshold openings, hence the differ- 
ences between the various methods to extract asymptotic 
observables from an ionization wave packet are particu- 
larly transparent in this case. 

We carried out two separate simulations, one starting 
from the ground state l-'^S (Is^) and one starting from 
the 2^S (ls2s) metastable state (lifetime 19.7 ms [88]) of 
the atom. A short (Tfwhm=200 as), moderately intense 
(Ipeak = 10^^ W/cm^) Gaussian XUV pulse with carrier 
frequency lo = 2.4 a. u. and uj = 1.65 a.u., respectively, 
was employed. In both simulations the energy of the ex- 
cited component at the end of the pulse was centered on 
the N — 2 threshold in order to populate the doubly 
excited states below the threshold. While the oscillator 
strength for exciting the doubly excited resonances 
from the ground state is weak, the excitation probability 
is strongly enhanced starting from the 2^S (ls2s) singly 
excited state. 
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E(au) 

FIG. 2. Total photoelectron spectrum of the wave packet cre- 
ated by the action of a sub-femtosecond pulse on the ground 
state of the helium atom and computed by projecting the wave 
packet on products of bound (Z=2) and continuum (Z=l) 
Coulomb functions. Since such products are not eigenstates 
of the field-free Hamiltonian, the spectrum is only approxi- 
mated and changes with time after the XUV pulse. As the 
doubly excited states populated by the pulse decay, charac- 
teristic Fano profiles build up. See text for more details. 



A. Projection on Coulomb functions vs projection 
on scattering states 

In Fig. 2 the square modulus |(lsi?p|5'(i)) [^ of the pro- 
jection on the product of a bound {Z — 2) and a con- 
tinuum {Z — 1) Coulomb function of the wave packet 
starting from the ground state is shown as a function of 
time. The photoelectron spectrum increases rapidly from 
zero to a smooth Gaussian profile immediately after the 
external field is over. At this stage, only the "direct ion- 
ization" component is visible, since the duration of the 
external pulse is much shorter than the lifetime of any of 
the metastable states, hence the smoothness of the 
profile right at the end of the pulse. As time passes, the 
localized part of the metastable states that were popu- 
lated by the pulse progressively decay and Fano-like in- 
terference features emerge as a consequence. The most 
prominent among these resonances, seen close to the cen- 
ter of the spectrum, are those converging to the N=2 
threshold and correspond mostly to the sp'^ series of the 
autoionizing states [89]. Three major features of this 
time-dependent spectrum are to be noted. First, each 
resonance peak takes a long time to reach convergence in 
comparison to the duration of the laser pulse. Second, 
decaying resonances give rise to quantum beats due to the 
interference between the direct and the delayed electron 
emission [90-92] . Such beats in the time domain are ob- 
servable when the wave packet is probed at finite times, 
e.g. by a pump-probe protocol. Third, a whole Rydberg 
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series of autoionizing states, like the one visible, com- 
prises states with arbitrarily long lifetimes which, with 
this method, cannot realistically ever be spectrally re- 
solved. A close-up of the spectrum near the spj ^P" (or 
2s2p) resonance (Fig. 3) shows the evolution of the pho- 
toelectron spectrum as a function of time determined by 
projection on the product of Coulomb functions and the 
asymptotic (t — >■ oo) limit of the Fano profile obtained by 
projecting the wave packet immediately after the end of 
the pulse on the tpi^^p scattering states. The scattering 
states account automatically for the time evolution of the 
full wave packet to infinity and no free propagation be- 
yond the end of the pulse is required. In turn, quantum 
beats occurring at finite time are not visible in the PSS 
method as it projects into the asymptotic future. When 
the projection is made on scattering states, the propaga- 
tion can be carried out in a box as small as 200 Bohr radii 
sufficient to contain the wave packet close to the end of 
the pulse. This is to be contrasted with the spectrum 
obtained by the projection on the product of Coulomb 
functions that has still not reached the asymptotic limit 
even though the box is an order of magnitude larger in 
size. The spectrum obtained by projection onto partial- 
wave channel functions with a frozen core (method (ii) 
in section IV) yields results very close to the projection 
onto Coulomb functions (Fig. 3) and will, therefore, not 
be separately discussed in the remainder of the paper. 



B. Conversion of B-spline close-coupling functions 
to the FEDVR basis 

Table I compares the energies of ^S*^, ^P°, and ^D*^ 
helium Rydberg states with principal quantum number 
n for the outer electron up to n — 6, obtained by diago- 
nalizing the full configuration interaction Hamiltonian of 
helium built in either the FEDVR or the B-spline basis. 
The FEDVR basis comprised eleven functions per finite 
element. The width of the first element was 2a.u. and 
increased linearly to 4.0 a. u. within the first 5 finite el- 
ements. The grid extensions were 28 a.u. for one radial 
coordinate and 156 a.u. for the other. The B-spline basis 
comprised spline functions of order k — 8 [93] defined 
on a non-uniform grid of nodes, optimized at small radii 
to optimize the representation of the ground state of the 
helium atom, and with an asymptotic spacing between 
consecutive nodes of 0.5 a.u. up to a maximum radius 
of 800 a.u. . For both the FEDVR and the B-spline ba- 
sis, the maximum orbital angular momentum Imax = 4 
(which suffices for the present comparison between meth- 
ods) was used. The very good agreement between the 
two approaches for the Rydberg spectrum indicates that 
the wave functions in the two bases are represented at 
comparable levels of accuracy. 

To assess the accuracy with which the wave functions 
computed in the B-spline basis are converted to the 
FEDVR basis we computed the norm of the converted 
Rydberg states {(l)n\4>n) as well as their overlap {(f>n\4>n) 




FIG. 3. Close-up of the evolution of the photoelectron 
spectrum determined by projection on products of Coulomb 
functions in comparison with the asymptotic Fano profile ex- 
tracted from the same wave packet by projection onto scat- 
tering states immediately after the end of the external pulse. 



with the corresponding states computed directly in the 
FEDVR basis. Both numbers should within the numer- 
ical accuracy be close to 1. The errors 5„ = 1 — 
and Sn = 1 — {(f>n\'Pn) for the states listed in table I are 
between 10"" and lO^'^ (table II). 

This confirms the accuracy of the conversion from the 
B-spline to the FEDVR basis. The error is larger for the 
ground state than for the excited states because the first 
finite element is still comparatively wide. This minor 
discrepancy could be easily fixed by choosing a smaller 
radial span for the first few finite elements. 

C. Comparison between methods 

The photoelectron spectrum from the ground state 
(Fig. 4) is computed (a) by projecting the wave func- 
tion on products of bound (Z=2) and continuum (Z=l) 
Coulomb functions after a lapse of time of 1590 as from 
the peak of the XUV pulse, (b) by projecting the wave 
function on the true scattering states of the field-free 
Hamiltonian, and (c) by extracting it with the Bcrkeley- 
ECS method. While the projection on Coulomb func- 
tions represents rather well the smooth background spec- 
trum associated to the direct ionization component, it 
fails to reproduce the sharp features associated to the 
autoionizing states. The spectra calculated with the PSS 
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TABLE I. Comparison between the energies of the He Ry- 
dberg states with principal quantum number for the outer 
electron up to n = 6, obtained by diagonalizing the Hamilto- 
nian in either the FEDVR (upper value, R — 156 a.u.) or the 
B-spline (lower value, i?=800 a.u.) basis. 





Symmetry 






n ^S" 


Ipo 






1 -2.903 5102 








-2.903 5164 








2 -2.145 9610 


-2.123 8231 






-2.145 9615 


-2.123 8232 






3 -2.061 2684 


-2.055 1399 




-2.055 6203 


-2.061 2685 


-2.055 1400 




-2.055 6203 


4 -2.033 5852 


-2.031 0669 




-2.031 2796 


-2.033 5853 


-2.031 0669 




-2.031 2796 


5 -2.021 1761 


-2.019 9046 




-2.020 0016 


-2.021 1761 


-2.019 9045 




-2.020 0016 


6 -2.014 5627 


-2.013 8331 




-2.013 8981 


-2.014 5627 


-2.013 8331 




-2.013 8981 


TABLE II. Error in the norm of the bound states translated 


from the B-spline to the FEDVR basis Sn = 


1- 


((^„ (^„), and 


error in the overlap 5„ 


= 1 - {</>n (^„) (see 


text for details). 


The notation [n] is a shorthand for 10 








Ipo 






n Sn Sn 


^71 Sn 


Sn Sn 


1 8.0[-7] 8.2[-7] 








2 7.6[-8] 7.8[-8] 


7.4[-ll] 3.9[-10] 






3 2.4[-8] 2.5[-8] 


2.5[-9] 2.7[-9] 


3.3| 


-9] 3.3[-9] 


4 1.2 [-8] 2.1 [-8] 


3.0[-9] 1.3[-8] 


3.3| 


-9] 9.8[-9] 


5 7.6[-9] 1.7[-8] 


3.1[-9] l.l[-8] 


3.3[-9] 1.3[-8] 


6 6.1[-9] 1.3[-8] 


4.7[-9] 1.6[-8] 


4.0[-9] 1.3[-8] 



and with the Berkeley-ECS method feature accurately 
a large number of resonant profiles and are in excel- 
lent agreement with each other. We emphasize that the 
smaller number of resonances appearing in (c) than in 
(b) is not due to a fundamental limitation of the method 
but because of the use of a coarser energy grid. Addi- 
tional insights can be gained from a close-up (Fig. 5a) and 
a logarithmic presentation of the photoionization prob- 
ability (Fig. 5b). Fig. 5 highlights two aspects of the 
performance of the three methods. First, the projection 
on scattering states and the Berkeley-ECS method are 
in excellent agreement close to the resonance (Fig. 5a) 
and down to the smallest probability densities (Fig. 5b). 
The deviation of the background profile of the spectrum 
from a parabola at low energies is due to the fact that 
the Gaussian envelope of the XUV pulse is eventually 
truncated. Second, the spectrum obtained through the 
projection on the product of Coulomb functions clearly 
deviates from the background below —1 a.u. , however, 
only when P{E) is already small (^ 10~^ of the direct 
ionization peak). This error due to the contamination by 
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FIG. 4. Photoelectron distribution in the Is channel 
resulting from the photoionization of helium from the ground 
state with an XUV pulse with uj = 2.4 a.u., / = 10^^W/cm^ 
and rFWHM=200 as. Projection on scattering states (PSS), 
red solid line; Berkeley-ECS, blue dashed line; projection on 
the product of Coulomb functions a.t t = 1590 as after the 
center of the XUV pulse, green dotted line. 



doubly excited states is small in the present case since 
they provide only a minor admixture to the wave packet. 
When the metastable states have a higher relative weight, 
however, their spurious effect on the spectrum will be 
more significant. 

The convergence of the DES spectrum as a function of 
the size of the close-coupling expansion within the PSS 
is illustrated in Fig. 6. Here, the PSS employs a mini- 
mal close-coupling expansion involving the open N — 2 
channels (IsCp, 2sep, 2pes, and 2ped) only. Clearly, since 
the N=3 channels are not included, the higher mem- 
bers of the autoionizing Rydberg series converging to the 
N=3 threshold are not reproduced. As a consequence, 
the spectrum obtained with the PSS deviates from the 
one obtained with the Berkeley-ECS method at energies 
higher than -0.3 a.u. . On the other hand, the spectrum 
below the = 3 threshold is already well converged, 
i.e. the influence of the closed channels are adequately ac- 
counted for. This observation highlights the salient fea- 
ture of the close-coupling expansion with pseudostates, 
namely to drastically truncate the representation while 
still obtaining an accurate representation of the contin- 
uum states in any given single-ionization energy region. 

The emission spectrum from the metastable ls2s ^S 
excited state of helium ionized by an XUV pulse with 
central frequency uj = 1.65 a.u. (Fig. 7) reflects the much 
enhanced role doubly excited states play. Consequently, 
the errors resulting from projection onto Coulomb func- 
tions are more prominent (Fig. 7) . Even though the back- 
ground contribution is semi-quantitatively reproduced, 
the resonant part which dominates the spectrum is not. 
Moreover, at low energies the background, if small, is 
strongly overestimated and as much as two orders of 
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Photoelectron distribution 
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FIG. 5. Photoelectron spectrum as in Fig. 4 with (a) a 
close-up near the spj resonance, and (b) the spectrum on a 
logarithmic scale to highlight small deviations in the tails of 
the direct ionization peak. 



magnitude larger in absolute terms than in the case of 
the photoionization from the ground state. The spec- 
tra from the other two methods show again an excellent 
agreement. To complete the comparison between the 
Berkeley-ECS method and the PSS, we verify whether 
the two methods predict also the same photoelectron an- 
gular distributions. In contrast to energy distributions, 
angular distributions depend on the relative phases be- 
tween photoionization amplitudes in different channels 
as well as on their absolute values. It is thus sufficient 
to compare the absolute phases of the photoionization 
amplitudes in the different channels to ascertain whether 
the angular distribution would also agree. As an example 
of such a comparison, we present the partial photoelec- 
tron distribution in the 2pes channel from the ls2s state 
(Fig. 8a), and the asymptotic phase of the correspond- 
ing amplitude in the interaction representation (Fig. 8b) . 
We show here the PSS with two different truncations of 
the close-coupling expansion, including and excluding the 
= 3 channels. The agreement in both probability den- 
sity and phase between the Berkeley-ECS and the PSS is 
excellent for all open channels included in the expansion. 
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FIG. 6. Photoelectron distribution in the N — 2 ^P° channel 
resulting from the photoionization of helium from the ground 
state with an XUV pulse with uj = 2.4 a.u., / = lO'^W/cm^, 
and Tfwhm=200 as. PSS: red dashed line; ECS: blue solid 
line. 
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FIG. 7. Photoelectron distribution in the Is channel, fol- 
lowing photoionization of the 2^S metastable state of he- 
lium by an XUV pulse with central frequency uj — 1.65 a.u., 
/ = lO^^W/cm^, and rpwHM=200 as. PSS: red solid line; 
Berkeley-ECS: blue dotted line; projection onto Coulomb 
waves: green dashed line. 



D. Double ionization spectra 

To complete our comparison for the applicability and 
validity of alternative methods we analyze the continuum 
components of a correlated wave packet for the doubly 
ionized part of the spectrum caused by two-photon ab- 
sorption of the previously used Gaussian XUV pulse with 
7FWHM=200as and a carrier frequency w = 2.4 a.u. from 
the helium ground state. In this regime, the PSS is 
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-0.5 -0.4 -0.3 -0.2 

Energy (au) 

FIG. 8. Partial photoelectron distribution and phase in the 
2pes channel from photoionization of the ls2s^S state of he- 
lium by an XUV pulse with lu — 1.65 a.u., / = lO^^W/cm^, 
and Tfwhm=200 as. Three different extraction methods are 
compared: PSS with a close coupling basis which includes 
the N=2 (dark yellow solid line) and N=3 channels (red thick 
solid line); Berkeley-ECS method (blue dashed line), (a) pho- 
toemission probability density P, (b) asymptotic phase of the 
photoelectron amplitude. 



not applicable because of the lack of accurate scattering 
states for the double continuum. However, the Berkeley- 
ECS method is able to impose the correct boundary con- 
ditions and obtain the asymptotic spectral information 
of a doubly ionized wave packet directly after the com- 
pletion of the laser pulse. The projection of the wave 
packet on an uncorrelated, symmetrized product of two 
Coulomb functions with Z — 2 is straightforward. How- 
ever, it requires the propagation of the wave packet to 
large distances in order to control and minimize the er- 
ror due to the neglect of the electron-electron interac- 
tion. In practice this yields accurate results as long as 
the box and angular momentum basis are large enough 
to correctly represent the wave function at the time of 
projection. We compare the singly differential photo- 
electron spectrum with the prominent two peaks of a 
sequential two-photon double ionization process. The 
projection onto Coulomb waves 7500 as after the peak 
of the pulse shows excellent agreement with the spec- 
trum obtained by the Berkeley-ECS method using the 
wave packet immediately after the conclusion of the pulse 




Energy (au) 



FIG. 9. Singly differential photoelectron distribution 

after two-photon double ionization of helium by a sub- 
femtosecond XUV pulse (a; = 2.4 a.u., / = lO^^W/cm^, and 
rFWHM=200 as). Two different extraction methods are com- 
pared: projection on products of Coulomb functions with 
Z = 2 performed 15Q0as (red solid line) and 7500as (black 
solid line) after the peak of the XUV pulse and the Berkeley- 
ECS method (green dashed line) for a computational box with 
R = 244 a.u. performed ISOOas after the peak of the XUV 
pulse. 



(Fig. 9). In fact, virtually the same level of agreement 
is already reached for considerably smaller propagation 
times before the projection (not shown). Even projection 
directly after the conclusion of the field (1500 as after the 
peak of the XUV pulse, red line in Fig. 9) gives very 
similar results. We also find very good agreement for 
the angular distributions, e.g. for the energy-integrated 
conditional angular emission probability for one electron 
when the other electron is ejected along the laser polar- 
ization axis (Fig. 10a). The remaining small residual dif- 
ferences between the two methods can be understood by 
inspecting their respective convergence behaviour: The 
projection on Coulomb functions becomes more accurate 
when the electrons are propagated further apart before 
the spectral analysis is performed, provided the partial 
wave expansion of the wave function covers enough angu- 
lar momenta to accurately describe the motion of the free 
electrons. For the present case convergence is reached for 
propagation times of about 7500 as after the peak of the 
XUV pulse (see the solid hnes in Fig. 10b). The ex- 
tension of the computational box for the transformation 
to the spectral domain (i.e. the solution of Eq. 15) and 
for the surface integral (cf. Eq. 18) influences the quality 
of the results. Thus, for accurate angular distributions 
comparably large radial boxes are also required for the 
Berkeley-ECS method (see the dashed lines in Fig. 10b). 
The convergence behaviour of the two methods is linked: 
their results for the angular distributions agree when the 
extraction radius of the Berkeley-ECS method roughly 
equals the position of the (relevant parts) of the wave 
packet at the time of projection on the Coulomb func- 
tions (compare solid and dashed lines in Fig. 10b). Thus, 
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for larger ECS boxes even the small differences to the 
(converged) Coulomb projection in Fig. 10a would van- 
ish. For both methods the convergence for the angular 
distributions is considerably slower than for the energy 
spectra, especially where both electrons are emitted in 
the same direction (Fig. 10b). 
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FIG. 10. Conditional angular distribution for one electron 
emitted in the laser polarization axis [grey arrow in (a)] and 
integrated over both electron energies after two-photon dou- 
ble ionization of helium by a sub-femtosecond XUV pulse 
(u = 2.4a.u., J = lO^^W/cm^, and rFWHM=200 as). Two 
different extraction methods are compared (a): projection on 
products of Coulomb functions with Z = 2 performed 7500 as 
(solid black line) after the peak of the XUV pulse and the 
Berkeley-ECS method (dashed green line) for an extraction 
radius of J? = 244 a.u. performed 1500 as after the peak of the 
XUV pulse. In (b) a close-up of the emission of the two elec- 
trons in the same direction is shown. The different lines illus- 
trate the convergence of the Berkeley-ECS method as a func- 
tion of the extraction radius (dashed lines) and the Coulomb 
projection as a function of the propagation time (solid lines). 
The two techniques agree when the extraction radius of the 
Berkeley-ECS method approximately equals the position of 
the (relevant parts) of the wave packet at the time of projec- 
tion on the Coulomb functions, see text. In both cases, the 
partial wave expansion of the wave function includes angular 
momenta up to h = h = 15. 



VI. CONCLUSIONS 

We have analyzed and evaluated different methods to 
extract the continuum component of a correlated multi- 
electron wave function in helium obtained from the solu- 
tion of the time-dependent Schrodinger equation. Each 
of the three methods investigated, the projection onto 
Coulomb continuum states, the projection onto exact 
scattering states (PSS), and the Berkeley-ECS method 
feature both advantages and disadvantages depending on 
the energy range of interest, the required box size, and 
whether single or double ionization distributions are de- 
sired. Below the double ionization threshold, methods 
to extract photoelectron spectra that are based on the 
projection on approximated continuum channels such as 
Coulomb continuum functions work well when no res- 
onances are significantly populated, but become inap- 
plicable when spectral regions in the vicinity of narrow 
resonances, are populated. These methods converge ex- 
tremely slowly, in particular close to thresholds when 
the wave packet requires long times to reach the asymp- 
totic region. In contrast, remarkably good agreement 
between the PSS and the Berkeley-ECS method is ob- 
served. This finding suggests an efficient way of moni- 
toring the composition of a complex wave packet ip in 
the single-ionization channels below the double ioniza- 
tion threshold. Since the scattering states can be com- 
puted separately, once and for all, in an adapted basis, 
the extraction of the expansion coefficients of ip requires 
only the calculation of a simple scalar product. More- 
over, in contrast to the method based on exterior com- 
plex scaling, it opens the way to isolate portions of the 
wave packet by constructing projections onto specific pre- 
selected channels. Our current implementation exploits 
the fact that correlated single-ionization multi-channel 
scattering states of helium computed on a B-spline ba- 
sis can be accurately translated to a FEDVR basis in 
which the wave packet propagation is performed. Such 
an approach could be extended to alternative methods 
to accurately compute multi-channel scattering states for 
atoms and molecules (for example the R-matrix method) 
in order to extract the wave packet information following 
the integration of the TDSE on a grid. Above the dou- 
ble ionization threshold, the PSS becomes inapplicable 
while the Berkeley-ECS method agrees well with projec- 
tion onto Coulomb functions provided the computational 
box is sufficiently large and the wave packet is propagated 
into the asymptotic region. 
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